Analiza Nutrycji w Menu Starbucks: Badanie Rozkładów, Korelacji i Klasyfikacji

Opis danych i cel projektu

Dane są informacjami o wartościach odżywczych pozycji w menu sieci kawiarni Starbucks. Napoje podzielone są na kategorie:
(1) Kawa,
(2) Klasyczne napoje espresso,
(3) Charakterystyczne napoje espresso,
(4) Napoje herbaciane Tazo®,
(5) Mrożone napoje wstrząśnięte,
(6) Smoothies,
(7) Kawa mieszana Frappuccino®,
(8) Kawa mieszana Frappuccino® Light,
(9) Frappuccino® Blended Crème.

Zmienne V1, V2, … , V15 określają kolejno zawartość: kalorii, tłuszczu ogółem (g), tłuszczu trans (g), tłuszczu nasyconego (g), sodu (mg), węglowodanów ogółem (g), cholesterolu (mg), błonniku pokarmowego (g), cukrów (g), białka (g), witaminy A (% DV), witaminy C (% DV), wapnia (% DV), żelaza (% DV) i kofeiny (mg).

Źródło danych: https://www.kaggle.com/datasets/henryshan/starbucks.

Celem projektu jest klasyfikacja bez nadzoru i pod nadzorem danych na podstawie wyselekcjonowanych 15 wartości odżywczych.

Plan pracy:

1. Określenie podstawowych parametrów rozkładów 15 zmiennych
   
2. Zbadanie obserwacji odstajacych i wielomodalnosci rozkładów 15 zmiennych
   
3. Macierz korelacji, redukcja wymiaru przy użyciu pca i t-sne
   
4. Analiza skupień
   
5. Klasyfikacja z użyciem lasów losowych
   
6. Klasyfikacja z użyciem knn
   
7. Klasyfikacja z użyciem naiwnego klasyfikatora Bayes’a
   
8. Wnioski końcowe

Analiza danych jednowymiarowych

starbucks <- read.csv('starbucks_data.csv')
head(starbucks)

##   Category  V1  V2  V3  V4 V5 V6 V7 V8 V9 V10 V11 V12  V13 V14 V15
## 1        1   3 0.1 0.0 0.0  0  5  0  0  0 0.3 0.0   0 0.00   0 175
## 2        1   4 0.1 0.0 0.0  0 10  0  0  0 0.5 0.0   0 0.00   0 260
## 3        1   5 0.1 0.0 0.0  0 10  0  0  0 1.0 0.0   0 0.00   0 330
## 4        1   5 0.1 0.0 0.0  0 10  0  0  0 1.0 0.0   0 0.02   0 410
## 5        2  70 0.1 0.1 0.0  5 75 10  0  9 6.0 0.1   0 0.20   0  75
## 6        2 100 3.5 2.0 0.1 15 85 10  0  9 6.0 0.1   0 0.20   0  75

any(is.na(starbucks))

## [1] FALSE

Pierwszym krokiem będzie prezentacja histogramów oraz podstawowych statystyk opisowych dla każdej zmiennej.

Podsumowanie:

apply(starbucks[, -1], 2, summary)

##               V1        V2       V3         V4        V5       V6       V7
## Min.      3.0000  0.000000 0.000000 0.00000000  0.000000   0.0000  0.00000
## 1st Qu. 122.5000  0.200000 0.100000 0.00000000  0.000000  80.0000 21.00000
## Median  190.0000  2.500000 0.500000 0.00000000  5.000000 127.5000 36.00000
## Mean    201.0413  3.074312 1.399541 0.03944954  6.674312 136.8578 37.16055
## 3rd Qu. 270.0000  5.000000 2.000000 0.10000000 10.000000 180.0000 52.50000
## Max.    510.0000 15.000000 9.000000 0.30000000 40.000000 340.0000 90.00000
##                V8       V9       V10       V11        V12       V13        V14
## Min.    0.0000000  0.00000  0.000000 0.0000000 0.00000000 0.0000000 0.00000000
## 1st Qu. 0.0000000 19.00000  4.000000 0.0400000 0.00000000 0.1000000 0.00000000
## Median  0.0000000 33.00000  6.000000 0.0800000 0.00000000 0.2000000 0.03000000
## Mean    0.8440367 33.94037  7.305046 0.1031193 0.03738532 0.2168349 0.08041284
## 3rd Qu. 1.0000000 44.00000 10.000000 0.1500000 0.00000000 0.3000000 0.10000000
## Max.    8.0000000 84.00000 20.000000 0.5000000 1.00000000 0.6000000 0.50000000
##               V15
## Min.      0.00000
## 1st Qu.  51.25000
## Median   75.00000
## Mean     89.93119
## 3rd Qu. 143.75000
## Max.    410.00000

Odchylenie standardowe:

sort(apply(starbucks[, -1], 2, sd), decreasing = TRUE)

##           V1           V6          V15           V7           V9           V5 
## 102.27230816  80.42698815  64.58923623  20.87723395  19.91320160   8.84163698 
##          V10           V2           V3           V8          V12          V13 
##   4.79903238   3.00748387   1.68767551   1.44435537   0.15070662   0.14558188 
##          V14          V11           V4 
##   0.10797386   0.08214517   0.07313689

Wariancja:

sort(apply(starbucks[, -1], 2, var), decreasing = TRUE)

##           V1           V6          V15           V7           V9           V5 
## 1.045963e+04 6.468500e+03 4.171769e+03 4.358589e+02 3.965356e+02 7.817454e+01 
##          V10           V2           V3           V8          V12          V13 
## 2.303071e+01 9.044959e+00 2.848249e+00 2.086162e+00 2.271249e-02 2.119408e-02 
##          V14          V11           V4 
## 1.165835e-02 6.747829e-03 5.349004e-03

Skośność:

sort(apply(starbucks[, -1], 2, skewness), decreasing = TRUE)

##       V12        V8       V11        V4        V5        V3       V14        V2 
## 5.3854396 2.7564703 1.8490660 1.7734325 1.6086669 1.5706864 1.5425005 1.0636256 
##       V15       V13       V10        V6        V9        V1        V7 
## 0.8687843 0.6600035 0.6542657 0.4530970 0.4330222 0.3647753 0.3533467

Wniosek: Największymi wariancjami cechują się zmienne V1, V6 oraz V15, które jednocześnie mają największe średnie. Wszystkie zmienne mają dodatnią skośność.

Histogramy:

for(i in 1:15) hist(starbucks[,i+1], main=colnames(starbucks)[i+1], xlab="wartość", ylab="częstość")

Wniosek: Histogramy wyraźnie pokazują dodatnią skośność zmiennych. Rozkład normalny przypominają jedynie zmienne V1 i V7 - wykonam dla nich test Shapiro-Wilka, a dla zmiennych V3-4 i V8-12 wykonam test na jednomodalność.

shapiro.test(starbucks$V1)

## 
##  Shapiro-Wilk normality test
## 
## data:  starbucks$V1
## W = 0.98489, p-value = 0.02008

shapiro.test(starbucks$V7)

## 
##  Shapiro-Wilk normality test
## 
## data:  starbucks$V7
## W = 0.97972, p-value = 0.003147

Wniosek: Żadna ze zmiennych nie może mieć rozkładu normalnego.

silverman.test(starbucks$V3, k = 1)

## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0.2872873

silverman.test(starbucks$V4, k = 1)

## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0.008008008

silverman.test(starbucks$V8, k = 1)

## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0.002002002

silverman.test(starbucks$V9, k = 1)

## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0.5795796

silverman.test(starbucks$V10, k = 1)

## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0.6706707

silverman.test(starbucks$V11, k = 1)

## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0

silverman.test(starbucks$V12, k = 1)

## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0

Wniosek: Przyjmując kryterium 5% błędu otrzymujemy informację, że zmienne V4, V8, V11 i V12 są wielomodalne. Zobaczmy jak wyglądają wykresy wartości p dla tych zmiennych:

silverman.plot(starbucks$V4)

silverman.plot(starbucks$V8)

silverman.plot(starbucks$V11)

silverman.plot(starbucks$V12)

Wniosek: Z wykresów wynika (nadal przyjmując kryterium 5% błędu), że zmienna V4 powinna być czteromodalna, V8 dwumodalna, V11 dwumodalna a V12 trójmodalna lub czteromodalna.

Kolejny test - Grubbsa - określi czy próba posiada wartości odstające.

apply(starbucks[, -1], 2, function(x) grubbs.test(x))

## $V1
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 3.02094, U = 0.95775, p-value = 0.2494
## alternative hypothesis: highest value 510 is an outlier
## 
## 
## $V2
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 3.96534, U = 0.92721, p-value = 0.005931
## alternative hypothesis: highest value 15 is an outlier
## 
## 
## $V3
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 4.50351, U = 0.90611, p-value = 0.0004397
## alternative hypothesis: highest value 9 is an outlier
## 
## 
## $V4
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 3.56250, U = 0.94124, p-value = 0.03307
## alternative hypothesis: highest value 0.3 is an outlier
## 
## 
## $V5
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 3.76918, U = 0.93423, p-value = 0.01402
## alternative hypothesis: highest value 40 is an outlier
## 
## 
## $V6
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 2.52580, U = 0.97047, p-value = 1
## alternative hypothesis: highest value 340 is an outlier
## 
## 
## $V7
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 2.53096, U = 0.97034, p-value = 1
## alternative hypothesis: highest value 90 is an outlier
## 
## 
## $V8
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 4.95443, U = 0.88636, p-value = 3.728e-05
## alternative hypothesis: highest value 8 is an outlier
## 
## 
## $V9
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 2.51389, U = 0.97074, p-value = 1
## alternative hypothesis: highest value 84 is an outlier
## 
## 
## $V10
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 2.6453, U = 0.9676, p-value = 0.8416
## alternative hypothesis: highest value 20 is an outlier
## 
## 
## $V11
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 4.83146, U = 0.89193, p-value = 7.508e-05
## alternative hypothesis: highest value 0.5 is an outlier
## 
## 
## $V12
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 6.38734, U = 0.81112, p-value = 2.027e-09
## alternative hypothesis: highest value 1 is an outlier
## 
## 
## $V13
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 2.63196, U = 0.96793, p-value = 0.8765
## alternative hypothesis: highest value 0.6 is an outlier
## 
## 
## $V14
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 3.88601, U = 0.93009, p-value = 0.008446
## alternative hypothesis: highest value 0.5 is an outlier
## 
## 
## $V15
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 4.95545, U = 0.88631, p-value = 3.706e-05
## alternative hypothesis: highest value 410 is an outlier

Wniosek: Przyjmując kryterium 5% błędu mamy, że jedynie zmienne V1, V6, V7, V9 i V10 nie mają wartości odstających.

Macierz korelacji i redukcja wymiaru

Dokonujemy obliczeń macierzy korelacji Pearsona oraz współczynnika zależności Kendalla, aby uniezależnić się od założenia normalności rozkładu danych.

heatmaply_cor(cor(starbucks[, -1]))

Wniosek: Mapa cieplna jest bardzo zróżnicowana. Najsilniej skorelowane są zbiory zmiennych {V1, V6, V7, V9}, {V4, V5}, {V10, V11, V13} oraz {V8, V12}. Słabiej skorelowane są zbiory takie jak {V3, V4, V5} i {V2, V3, V5, V10, V13}. Najsłabiej są skorelowane np. zmienne V7 i V14.

heatmaply_cor(cor(starbucks[, -1], method = 'kendall'))

Wniosek: Macierz korelacji Kendalla jest podobna do macierzy Pearsona, widać na niej jednak nieco słabsze skorelowanie między zmiennymi.\

prcomp(starbucks[, -1]) -> pca.starbucks
summary(pca.starbucks)

## Importance of components:
##                             PC1     PC2      PC3      PC4     PC5     PC6
## Standard deviation     126.1494 65.4197 40.69682 11.77600 6.70131 3.09895
## Proportion of Variance   0.7218  0.1941  0.07512  0.00629 0.00204 0.00044
## Cumulative Proportion    0.7218  0.9159  0.99102  0.99731 0.99935 0.99979
##                            PC7     PC8     PC9    PC10    PC11    PC12    PC13
## Standard deviation     1.89625 0.84376 0.46864 0.42542 0.07908 0.05196 0.03425
## Proportion of Variance 0.00016 0.00003 0.00001 0.00001 0.00000 0.00000 0.00000
## Cumulative Proportion  0.99995 0.99998 0.99999 1.00000 1.00000 1.00000 1.00000
##                           PC14    PC15
## Standard deviation     0.02352 0.01781
## Proportion of Variance 0.00000 0.00000
## Cumulative Proportion  1.00000 1.00000

prcomp(starbucks[, -1], scale. = TRUE) -> pca.starbucks.sc
summary(pca.starbucks.sc)

## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.4715 1.5835 1.4521 1.16984 1.07465 0.88856 0.54014
## Proportion of Variance 0.4072 0.1672 0.1406 0.09124 0.07699 0.05264 0.01945
## Cumulative Proportion  0.4072 0.5744 0.7150 0.80620 0.88319 0.93583 0.95528
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.51786 0.43590 0.33292 0.24555 0.15849 0.11681 0.04168
## Proportion of Variance 0.01788 0.01267 0.00739 0.00402 0.00167 0.00091 0.00012
## Cumulative Proportion  0.97315 0.98582 0.99321 0.99723 0.99890 0.99981 0.99993
##                           PC15
## Standard deviation     0.03235
## Proportion of Variance 0.00007
## Cumulative Proportion  1.00000

df.pca <- data.frame(pc1 = pca.starbucks$x[, 1], pc2 = pca.starbucks$x[, 2], pc3 = pca.starbucks$x[, 3], kl = as.factor(starbucks$Category))
plot(x = df.pca[, 1], y = df.pca[, 2], col = as.factor(starbucks[, 1]))

plot_ly(df.pca, x=~pc1, y=~pc2, z=~pc3, color=~kl, type='scatter3d', colors = viridis(9))

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

pca.starbucks$rotation[, 1] -> pca.rot1
sort(pca.rot1)

##           V15            V4           V12           V11           V14 
## -7.695161e-05  1.752291e-04  1.761441e-04  2.217498e-04  3.216063e-04 
##           V13            V8            V3            V2           V10 
##  5.175786e-04  3.402807e-03  8.268791e-03  1.350765e-02  1.877526e-02 
##            V5            V9            V7            V6            V1 
##  2.470073e-02  1.443410e-01  1.536808e-01  5.770364e-01  7.882585e-01

pca.starbucks$rotation[, 2] -> pca.rot2
sort(pca.rot2)

##           V15            V6            V5            V2            V3 
## -9.796140e-01 -1.628696e-01 -5.026304e-03 -3.768940e-03 -2.980319e-03 
##           V13            V4           V14           V11           V10 
## -2.306455e-04 -5.301965e-06  9.095648e-06  2.818636e-05  3.950922e-04 
##           V12            V8            V9            V7            V1 
##  6.530707e-04  4.909092e-03  3.095668e-02  3.964840e-02  1.059556e-01

pca.starbucks$rotation[, 3] -> pca.rot3
sort(pca.rot3)

##            V6            V4           V11           V14           V13 
## -0.7984239836  0.0002377182  0.0002803016  0.0004738972  0.0007056207 
##           V12            V3            V8            V2           V10 
##  0.0007911645  0.0073478188  0.0094140046  0.0211193312  0.0300239710 
##            V5            V9            V7           V15            V1 
##  0.0322583642  0.0425003221  0.0735806015  0.1974442763  0.5601627919

Wniosek: Dla niewyskalowanego PCA mamy, że już pierwsze dwie składowe zawierają 92% wariancji, a trzy składowe 99%, natomiast dla wyskalowanego PCA trzy składowe zawierają zaledwie 72%. Dalej dla niewyskalowanego PCA: PC1 i PC2 zdominowane są przez zmienną V15. Niestety na wykresach nie widać aby ta zmienna oddzielała którykolwiek klaster od reszty.

set.seed(10)
Rtsne(starbucks[, -1], theta = 0, dims = 2) -> starbucks.tsne2
summary(starbucks.tsne2$Y)

##        V1                V2          
##  Min.   :-8.3469   Min.   :-16.0396  
##  1st Qu.:-5.0850   1st Qu.: -6.5213  
##  Median : 0.3785   Median :  0.1484  
##  Mean   : 0.0000   Mean   :  0.0000  
##  3rd Qu.: 4.1747   3rd Qu.:  6.9102  
##  Max.   : 9.8827   Max.   : 14.1478

sd(starbucks.tsne2$Y)

## [1] 6.805308

tsne2.y.df <- as.data.frame(starbucks.tsne2$Y)
ggplot(tsne2.y.df, aes(x = V1, y = V2, col = starbucks$Category)) +
  geom_point() + labs(x = "t-SNE 1", y = "t-SNE 2", col = "Category") + theme_classic()

Rtsne(starbucks[, -1], theta = 0, dims = 3) -> starbucks.tsne3


summary(starbucks.tsne3$Y)

##        V1                V2                 V3          
##  Min.   :-17.132   Min.   :-25.8961   Min.   :-22.5534  
##  1st Qu.: -5.677   1st Qu.:-10.0151   1st Qu.: -8.3242  
##  Median : -0.433   Median : -0.9437   Median :  0.5223  
##  Mean   :  0.000   Mean   :  0.0000   Mean   :  0.0000  
##  3rd Qu.:  8.073   3rd Qu.: 10.4523   3rd Qu.:  8.4342  
##  Max.   : 13.631   Max.   : 26.3293   Max.   : 17.9610

sd(starbucks.tsne3$Y)

## [1] 11.15663

tsne3.y.df <- as.data.frame(starbucks.tsne3$Y)


plot_ly(tsne3.y.df, x = ~V1, y = ~V2, z = ~V3, color = starbucks$Category, type = 'scatter3d', colors = "Set3")

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Wniosek: T-SNE jedynie trochę lepiej radzi sobie z odróżnianiem kategorii, nadal jest to raczej chaotyczna mieszanka.

Analiza skupień

Sprawdzę jak działa k-średnich dla ilości klastrów równej ilości kategorii kawy.

starbucks.km <- kmeans(starbucks[, -1], centers = 9)

df.pca2 <- data.frame(PC1 = pca.starbucks$x[, 1], PC2 = pca.starbucks$x[, 2], PC3 = pca.starbucks$x[, 3], kl = as.factor(starbucks.km$cluster))
plot_ly(df.pca2, x = ~PC1, y = ~PC2, z = ~PC3, color = starbucks.km$cluster, type = 'scatter3d')

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Wniosek: Wizualnie k-średnich ładnie podzieliło wykres.

Porównam jednak w tabeli czy jest to jakkolwiek wierne faktycznemu podziałowi.

table(starbucks$Category, starbucks.km$cluster)

##    
##      1  2  3  4  5  6  7  8  9
##   1  0  0  3  1  0  0  0  0  0
##   2 17 12  2  8  7 12  0  0  0
##   3  7  1  0  0  7  4  9  6  6
##   4 13  0  0  0  8  3  6  3  3
##   5  1  0  1  7  3  0  0  0  0
##   6  0  0  0  0  0  0  8  0  0
##   7  0  1  0  0  7  0  1 14 13
##   8  1  7  0  0  0  0  0  3  1
##   9  0  0  0  0  0  0  9  3  0

Wniosek: Jak widać z tabeli k-średnich słabo poradziło sobie z klasteryzacją, np. zmienna V2 została w głównej mierze podzielona między kl1 a kl5, podobnie zmienna V7 pomiędzy kl2 i kl7.

Użyję metody łokcia do znalezienia optymalnej ilości klastrów.

wss <- rep(NA, 9)
for(i in 1:9) wss[i] <- kmeans(starbucks[, -1], centers = i + 1)$tot.withinss
plot(1:9, wss)

Wniosek: Według wykresu optymalna ilość klastrów to 5.

Zatem teraz sprawdzę k-średnich dla ilości klastrów równej 5.

starbucks.km5 <- kmeans(starbucks[, -1], centers = 5)

df.pca3 <- data.frame(PC1 = pca.starbucks$x[, 1], PC2 = pca.starbucks$x[, 2], PC3 = pca.starbucks$x[, 3], kl = as.factor(starbucks.km$cluster))
plot_ly(df.pca3, x = ~PC1, y = ~PC2, z = ~PC3, color = starbucks.km5$cluster,type = 'scatter3d')

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

table(starbucks$Category, starbucks.km5$cluster)

##    
##      1  2  3  4  5
##   1  0  0  0  0  4
##   2 29  0  2 21  6
##   3 11 11  8 10  0
##   4 17  8  4  7  0
##   5  8  0  0  0  4
##   6  0  8  0  0  0
##   7  0  1 16 19  0
##   8  3  0  1  8  0
##   9  0 10  2  0  0

Wniosek: Nadal jest to bardzo nieprecyzyjny podział.

plot(agnes(dist(starbucks[, -1], method = 'minkowski', p = 1), diss = T))

plot(diana(starbucks[, -1]))

Wniosek: Metoda aglomeracyjna i podziałowa dają względem siebie bardzo zbliżone rezultaty.

Klasyfikacja pod nadzorem - lasy losowe

rfcv(trainx = starbucks[, -1], trainy = as.factor(starbucks$Category)) -> rf.starbucks
rf.starbucks

## $n.var
## [1] 15  8  4  1
## 
## $error.cv
##        15         8         4         1 
## 0.1513761 0.1238532 0.1376147 0.4357798 
## 
## $predicted
## $predicted$`15`
##   [1] 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 2 2 2 2 3 2 2 3 2 2 3 3 2 3 2
##  [38] 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2
##  [75] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4
## [112] 4 3 4 4 4 4 4 4 4 4 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5
## [149] 5 5 6 6 6 6 6 6 6 6 9 7 8 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## [186] 7 7 7 7 7 7 7 7 7 8 8 8 8 8 7 8 8 8 8 7 7 9 9 9 9 9 9 9 9 7 7 7 7
## Levels: 1 2 3 4 5 6 7 8 9
## 
## $predicted$`8`
##   [1] 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 2 3 2 3 2 2 2
##  [38] 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 2 3 2 2 2 3
##  [75] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 3 4 4 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4
## [112] 4 3 4 4 4 4 4 4 4 4 3 4 4 4 4 4 4 4 4 4 4 4 4 3 4 3 4 5 5 5 5 5 5 5 5 5 5
## [149] 5 5 6 6 6 6 6 6 6 6 9 7 8 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## [186] 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 7 7 9 9 9 9 9 9 9 9 9 9 7 9
## Levels: 1 2 3 4 5 6 7 8 9
## 
## $predicted$`4`
##   [1] 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 2 3 2 2 2
##  [38] 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 3 2 3 3 2 3 2 3
##  [75] 3 3 2 7 3 7 3 3 2 3 3 3 4 3 4 3 3 4 3 3 3 3 3 3 4 3 3 3 4 4 4 4 4 4 4 4 4
## [112] 4 4 4 4 4 4 4 4 4 4 4 4 4 7 4 4 4 4 4 4 4 4 4 4 4 9 4 5 5 5 5 5 5 5 5 5 5
## [149] 5 5 6 6 6 6 6 6 6 6 9 7 7 7 7 7 7 7 7 4 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## [186] 3 3 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 3 9 9 7 9 9 9 9 9 9 9 7 9
## Levels: 1 2 3 4 5 6 7 8 9
## 
## $predicted$`1`
##   [1] 2 5 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 7 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [38] 2 2 2 2 2 5 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [75] 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 4 3 4 3 3 3 4 9 4 4 4 4 4 7 7 7 2 2 7
## [112] 8 8 8 3 4 3 4 4 4 4 4 4 7 7 7 4 9 9 9 9 4 4 9 9 9 9 9 8 8 2 5 5 5 5 5 5 5
## [149] 5 5 6 6 6 4 4 9 9 4 4 7 7 7 4 2 4 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## [186] 2 2 2 7 7 7 7 7 7 7 7 4 7 4 7 7 5 4 7 7 5 9 9 9 4 4 4 4 9 4 4 9 4
## Levels: 1 2 3 4 5 6 7 8 9

Wniosek: Błąd walidacji krzyżowej jest najniższy przy 8 cechach - wybierzemy te 8 najlepszych cech.

randomForest(starbucks[, -1], as.factor(starbucks$Category), importance = T) -> rf.starbucks.imp
varImpPlot(rf.starbucks.imp)

Wniosek: Najlepsze są cechy o numerach 15, 9, 7, 6, 10, 13, 1 i 11.

summary(as.factor(starbucks$Category))

##  1  2  3  4  5  6  7  8  9 
##  4 58 40 36 12  8 36 12 12

model <- rpart(as.factor(starbucks$Category) ~ ., method = "class", data=starbucks[, -1])
printcp(model)

## 
## Classification tree:
## rpart(formula = as.factor(starbucks$Category) ~ ., data = starbucks[, 
##     -1], method = "class")
## 
## Variables actually used in tree construction:
## [1] V12 V13 V15 V3  V6  V7  V8  V9 
## 
## Root node error: 160/218 = 0.73394
## 
## n= 218 
## 
##        CP nsplit rel error  xerror     xstd
## 1 0.11875      0   1.00000 1.00000 0.040778
## 2 0.08750      2   0.76250 0.77500 0.045701
## 3 0.07500      3   0.67500 0.73750 0.045983
## 4 0.06875      4   0.60000 0.73125 0.046016
## 5 0.05625      5   0.53125 0.64375 0.046070
## 6 0.05000      6   0.47500 0.61250 0.045905
## 7 0.03750      8   0.37500 0.60000 0.045811
## 8 0.02500      9   0.33750 0.50000 0.044477
## 9 0.01000     11   0.28750 0.47500 0.043975

fancyRpartPlot(model, )

table(predict(model, starbucks, type="class"), starbucks$Category)

##    
##      1  2  3  4  5  6  7  8  9
##   1  0  0  0  0  0  0  0  0  0
##   2  4 47  2  9  0  0  0  0  0
##   3  0  7 29  4  0  0  0  0  0
##   4  0  4  8 23  0  0  0  0  0
##   5  0  0  0  0 12  0  0  0  0
##   6  0  0  1  0  0  6  0  0  0
##   7  0  0  0  0  0  2 35  4  0
##   8  0  0  0  0  0  0  0  8  0
##   9  0  0  0  0  0  0  1  0 12

set.seed(10)
randomForest(starbucks[, c(16, 10, 8, 7, 11, 14, 2, 12)], as.factor(starbucks$Category)) -> rf.starbucks.sel
rf.starbucks.sel$confusion

##   1  2  3  4  5 6  7  8  9 class.error
## 1 3  1  0  0  0 0  0  0  0  0.25000000
## 2 0 55  3  0  0 0  0  0  0  0.05172414
## 3 0  8 27  4  1 0  0  0  0  0.32500000
## 4 0  0  5 31  0 0  0  0  0  0.13888889
## 5 0  0  0  0 12 0  0  0  0  0.00000000
## 6 0  0  0  0  0 8  0  0  0  0.00000000
## 7 0  0  0  0  0 0 34  1  1  0.05555556
## 8 0  0  0  0  0 0  2 10  0  0.16666667
## 9 0  0  0  0  0 0  0  0 12  0.00000000

1 - mean(rf.starbucks.sel$confusion[, 10])

## [1] 0.8902405

Wniosek: Metodą lasów losowych osiągnęliśmy średnią dokładność na poziomie 89%, co stanowi ulepszenie w porównaniu z drzewami decyzyjnymi. Największy błąd występuje w klasie V3 i wynosi 33%.

Klasyfikacja pod nadzorem - KNN i Naiwny klasyfikator Bayes’a

party <- createDataPartition(starbucks[, 1], p = 0.7, list = F)
starbucks.test <- starbucks[-party,]
starbucks.train <- starbucks[party,]

prop.table(table(starbucks$Category))

## 
##          1          2          3          4          5          6          7 
## 0.01834862 0.26605505 0.18348624 0.16513761 0.05504587 0.03669725 0.16513761 
##          8          9 
## 0.05504587 0.05504587

prop.table(table(starbucks.train$Category))

## 
##          1          2          3          4          5          6          7 
## 0.01935484 0.26451613 0.19354839 0.15483871 0.05806452 0.03225806 0.16774194 
##          8          9 
## 0.05806452 0.05161290

preparty <- preProcess(starbucks.train[, -1], method=c("center", "scale"))
starbucks.train.sc <- predict(preparty, starbucks.train)
starbucks.test.sc <- predict(preparty, starbucks.test)

apply(starbucks.train.sc[,-1],2,sd)

##  V1  V2  V3  V4  V5  V6  V7  V8  V9 V10 V11 V12 V13 V14 V15 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1

summary(starbucks.train.sc)

##     Category           V1                 V2                V3         
##  Min.   :1.000   Min.   :-1.95136   Min.   :-1.0067   Min.   :-0.8251  
##  1st Qu.:2.000   1st Qu.:-0.74160   1st Qu.:-0.9407   1st Qu.:-0.7665  
##  Median :4.000   Median :-0.09705   Median :-0.1816   Median :-0.5322  
##  Mean   :4.335   Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.:7.000   3rd Qu.: 0.69624   3rd Qu.: 0.5609   3rd Qu.: 0.3466  
##  Max.   :9.000   Max.   : 3.07610   Max.   : 3.9438   Max.   : 4.4474  
##        V4                V5                V6                V7          
##  Min.   :-0.5387   Min.   :-0.7773   Min.   :-1.6748   Min.   :-1.76838  
##  1st Qu.:-0.5387   1st Qu.:-0.7773   1st Qu.:-0.7316   1st Qu.:-0.71544  
##  Median :-0.5387   Median :-0.1867   Median :-0.1535   Median :-0.04539  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.: 0.9010   3rd Qu.: 0.4039   3rd Qu.: 0.5767   3rd Qu.: 0.64859  
##  Max.   : 3.7804   Max.   : 3.3569   Max.   : 2.4631   Max.   : 2.53909  
##        V8                V9                V10               V11         
##  Min.   :-0.5670   Min.   :-1.69731   Min.   :-1.5470   Min.   :-1.2910  
##  1st Qu.:-0.5670   1st Qu.:-0.69279   1st Qu.:-0.6898   1st Qu.:-0.7847  
##  Median :-0.5670   Median :-0.09008   Median :-0.2613   Median :-0.2785  
##  Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.1654   3rd Qu.: 0.51263   3rd Qu.: 0.5958   3rd Qu.: 0.6075  
##  Max.   : 4.5598   Max.   : 2.52166   Max.   : 2.5244   Max.   : 5.0375  
##       V12               V13               V14               V15         
##  Min.   :-0.2256   Min.   :-1.5027   Min.   :-0.7346   Min.   :-1.4288  
##  1st Qu.:-0.2256   1st Qu.:-0.8045   1st Qu.:-0.7346   1st Qu.:-0.6249  
##  Median :-0.2256   Median :-0.1063   Median :-0.5555   Median :-0.2230  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.:-0.2256   3rd Qu.: 0.5919   3rd Qu.: 0.1613   3rd Qu.: 0.9024  
##  Max.   : 6.6581   Max.   : 2.6864   Max.   : 3.7449   Max.   : 3.8766

Wniosek: Podzieliłam dane na część treningową (70%) i testową (30%), tak że ich procentowy podział klas pozostał bez większych zmian. Następnie wyskalowałam/wystandaryzowałam je tak, że w danych treningowych odchylenie standardowe wynosi 1, a średnia 0.

Wykorzystam ten podział danych do KNN (k najbliższych sąsiadów), policzę również optymalne k dzięki pętli for.

set.seed(10)

długość <- length(starbucks.train.sc$Category)
max.dokładność <- 0
max.k <- 0
for (i in c(2:długość)){
  prognoza <- knn(starbucks.train.sc[, -1],starbucks.test.sc[, -1],starbucks.train.sc$Category,k = i)
  dokładność <- mean(prognoza == starbucks.test.sc$Category)
if (max.dokładność < dokładność){
  max.k = i
  max.dokładność = dokładność
}
}
c(max.k, max.dokładność)

## [1] 3.0000000 0.8095238

knn(starbucks.train.sc[, -1],starbucks.test.sc[, -1], starbucks.train.sc$Category, k = max.k) -> starbucks.knn
conf_matrix.knn <- table(starbucks.knn, starbucks.test.sc$Category)
conf_matrix.knn

##              
## starbucks.knn  1  2  3  4  5  6  7  8  9
##             1  1  1  0  0  0  0  0  0  0
##             2  0 12  4  1  0  0  0  0  0
##             3  0  4  5  4  0  0  0  0  0
##             4  0  0  1  7  0  0  0  0  0
##             5  0  0  0  0  3  0  0  0  0
##             6  0  0  0  0  0  3  0  0  0
##             7  0  0  0  0  0  0 10  0  1
##             8  0  0  0  0  0  0  0  3  0
##             9  0  0  0  0  0  0  0  0  3

mean(starbucks.knn == starbucks.test.sc$Category)

## [1] 0.7460317

Wniosek: Najlepsze k dla moich danych to 2, a na macierzy widać że knn średnio sobie poradziło.

Wykonam teraz Naiwny klasyfikator Bayes’a, również przy użyciu stworzonych wcześniej danych treningowych i testowych.

model.nb <- naiveBayes(Category~., data = starbucks.train.sc)
model.nb

## 
## Naive Bayes Classifier for Discrete Predictors
## 
## Call:
## naiveBayes.default(x = X, y = Y, laplace = laplace)
## 
## A-priori probabilities:
## Y
##          1          2          3          4          5          6          7 
## 0.01935484 0.26451613 0.19354839 0.15483871 0.05806452 0.03225806 0.16774194 
##          8          9 
## 0.05806452 0.05161290 
## 
## Conditional probabilities:
##    V1
## Y         [,1]        [,2]
##   1 -1.9414457 0.009916107
##   2 -0.5771346 0.793167598
##   3  0.5078326 1.076330422
##   4 -0.1301036 0.739552706
##   5 -0.8242311 0.451700157
##   6  0.7953997 0.140234926
##   7  0.7763303 0.789819849
##   8 -0.2182468 0.577730830
##   9  0.3243846 0.692608034
## 
##    V2
## Y          [,1]      [,2]
##   1 -0.97370076 0.0000000
##   2  0.05181320 1.0076355
##   3  0.72046333 1.2669437
##   4 -0.12524358 0.7538134
##   5 -0.77568158 0.2640256
##   6 -0.34664003 0.4811001
##   7 -0.05341932 0.8134777
##   8 -0.51165601 0.6475695
##   9 -0.38789402 0.7400268
## 
##    V3
## Y         [,1]      [,2]
##   1 -0.8250782 0.0000000
##   2 -0.1449416 0.8320957
##   3  0.7664320 1.3069830
##   4 -0.2905065 0.6331589
##   5 -0.7014026 0.1974620
##   6 -0.4032792 0.3055326
##   7  0.1843554 1.0532350
##   8 -0.1806630 0.8884582
##   9 -0.3051523 0.6653312
## 
##    V4
## Y          [,1]      [,2]
##   1 -0.53873064 0.0000000
##   2  0.05822278 1.1140616
##   3  0.18112496 1.2396401
##   4  0.12113699 1.0381182
##   5 -0.37876273 0.4799037
##   6 -0.53873064 0.0000000
##   7  0.01500443 0.8222173
##   8 -0.53873064 0.0000000
##   9  0.18112496 1.0883194
## 
##    V5
## Y         [,1]      [,2]
##   1 -0.7773001 0.0000000
##   2  0.1157958 1.1215251
##   3  0.4038912 1.2789710
##   4  0.1332015 1.0732825
##   5 -0.5804349 0.4176142
##   6 -0.1867044 0.0000000
##   7 -0.2321349 0.6233082
##   8 -0.5148131 0.3112712
##   9 -0.1128800 0.6650058
## 
##    V6
## Y         [,1]       [,2]
##   1 -1.5733828 0.03513289
##   2 -0.3891940 0.62695087
##   3  0.0290519 0.97821058
##   4 -0.5388996 0.46157812
##   5 -1.3042820 0.23068517
##   6 -0.1048224 0.14528204
##   7  1.1337489 0.73016100
##   8  1.1649550 0.70441227
##   9  0.6299649 0.65574860
## 
##    V7
## Y          [,1]      [,2]
##   1 -1.76838012 0.0000000
##   2 -0.79949034 0.6654875
##   3  0.29442124 0.9239913
##   4 -0.15906003 0.6112599
##   5 -0.57185961 0.4068170
##   6  0.85439284 0.1490617
##   7  1.10142830 0.7963305
##   8  0.02374171 0.5340885
##   9  0.63662611 0.6767333
## 
##    V8
## Y          [,1]      [,2]
##   1 -0.56701657 0.0000000
##   2 -0.06684342 0.7220519
##   3 -0.02992587 0.6916995
##   4 -0.23133489 0.5705570
##   5 -0.56701657 0.0000000
##   6  4.26679970 0.4011500
##   7 -0.03180381 0.6059957
##   8 -0.07875230 0.6342739
##   9 -0.47546702 0.2589412
## 
##    V9
## Y          [,1]      [,2]
##   1 -1.69730775 0.0000000
##   2 -0.85694449 0.5996699
##   3  0.25982396 0.9215496
##   4 -0.08589691 0.6177917
##   5 -0.45282422 0.4222984
##   6  0.18113689 0.2058645
##   7  1.18874356 0.7714827
##   8  0.07175628 0.4763366
##   9  0.76375608 0.7072493
## 
##    V10
## Y          [,1]       [,2]
##   1 -1.41840270 0.07726015
##   2  0.32929493 0.91010271
##   3  0.59583973 1.13944092
##   4 -0.02914684 0.65030504
##   5 -1.08983833 0.30379590
##   6  2.01009505 0.28748828
##   7 -0.63215586 0.31816346
##   8 -0.57080186 0.30513643
##   9 -0.63627665 0.27465264
## 
##    V11
## Y          [,1]      [,2]
##   1 -1.29102775 0.0000000
##   2  0.32661582 0.8911882
##   3  0.39659023 0.8802654
##   4  0.08543566 0.6489549
##   5 -0.89725023 0.2861496
##   6  1.59479899 3.1444766
##   7 -0.59001721 0.2801336
##   8 -0.47534573 0.2109522
##   9 -0.46831399 0.2620276
## 
##    V12
## Y          [,1]       [,2]
##   1 -0.22560925 0.00000000
##   2 -0.21217755 0.04135841
##   3 -0.17053928 0.07754429
##   4 -0.05064735 0.36310884
##   5 -0.22560925 0.00000000
##   6  4.93720085 2.40931138
##   7 -0.22560925 0.00000000
##   8 -0.22560925 0.00000000
##   9  0.11857809 0.22077105
## 
##    V13
## Y         [,1]      [,2]
##   1 -1.5026628 0.0000000
##   2  0.4624575 0.9513657
##   3  0.7082394 1.1680382
##   4  0.1846047 0.7611138
##   5 -1.0372097 0.2792719
##   6 -0.6648472 0.3122354
##   7 -0.6594766 0.2792182
##   8 -0.6648472 0.2035526
##   9 -0.5426658 0.3094351
## 
##    V14
## Y          [,1]      [,2]
##   1 -0.73463724 0.0000000
##   2  0.05200585 1.1435463
##   3  0.19709780 1.1580197
##   4 -0.30535227 0.5751615
##   5 -0.61518403 0.2003291
##   6  0.55545742 0.9258587
##   7  0.24396021 1.0245962
##   8  0.20107957 1.1750555
##   9 -0.48826499 0.1900489
## 
##    V15
## Y         [,1]      [,2]
##   1  2.6708264 1.2478967
##   2  0.5553523 0.7936152
##   3 -0.2390779 0.9132074
##   4 -0.7655937 0.6418702
##   5  0.6880187 0.5131377
##   6 -1.3323015 0.1320845
##   7  0.1572776 0.5354125
##   8  0.2860982 0.5019985
##   9 -1.4287624 0.0000000

pred.nb <- predict(model.nb, starbucks.test.sc)
pred.nb

##  [1] 1 2 2 8 8 2 3 5 2 3 3 5 5 2 2 5 2 2 2 2 2 4 3 3 8 3 3 9 8 3 3 4 4 4 4 4 4 4
## [39] 5 8 1 1 5 6 6 6 8 7 7 8 8 8 8 7 7 7 8 8 8 9 9 9 8
## Levels: 1 2 3 4 5 6 7 8 9

conf_matrix.bayes <- table(pred.nb, starbucks.test.sc$Category)
conf_matrix.bayes

##        
## pred.nb 1 2 3 4 5 6 7 8 9
##       1 1 0 0 0 2 0 0 0 0
##       2 0 8 3 0 0 0 0 0 0
##       3 0 3 4 2 0 0 0 0 0
##       4 0 0 1 7 0 0 0 0 0
##       5 0 4 0 1 1 0 0 0 0
##       6 0 0 0 0 0 3 0 0 0
##       7 0 0 0 0 0 0 5 0 0
##       8 0 2 1 2 0 0 5 3 1
##       9 0 0 1 0 0 0 0 0 3

mean(pred.nb == starbucks.test.sc$Category)

## [1] 0.5555556

Wniosek: Naiwny klasyfikator Bayes’a również poradził sobie jedynie w około 60%.

Wnioski końcowe

Celem projektu było przeprowadzenie klasyfikacji danych zarówno bez nadzoru, jak i pod nadzorem, opierając się na 15 wartościach odżywczych.

Przeprowadziłam analizę rozkładów poszczególnych zmiennych, identyfikując wariancję, skośność i wielomodalność:
- wszystkie zmienne wykazują dodatnią skośność oraz nie są podobne do rozkładu normalnego,
- zmienne V4, V8, V11 i V12 są prawdopodobnie wielomodalne,
- 2/3 zmiennych posiada wartości odstające.

Zbadałam macierz korelacji Pearsona i współczynnika zależności Kendalla oraz zastosowałam PCA i t-SNE do redukcji wymiaru danych:
- macierze wykazały, że spora część zmiennych wykazuje silną korelację,
- metody redukcji wymiaru niestety nie zapewniły jednoznacznej separacji danych.

Przeprowadziłam analizę skupień za pomocą k-średnich, aglomeracyjnej i podziałowej, z identyfikacją optymalnej liczby klastrów. Wyniki analizy skupień były jednak niezadowalające, co sugeruje trudności w jednoznacznym podziale danych na grupy.

Zastosowałam algorytm lasów losowych do klasyfikacji, wykorzystując istotność cech do wyboru najlepszych predyktorów. Uzyskałam dokładność klasyfikacji na poziomie 89%.

Przeprowadziłam klasyfikację za pomocą k najbliższych sąsiadów (KNN) i naiwnego klasyfikatora Bayes’a. Oba modele osiągnęły dokładność na poziomie około 60%.

Pomimo zastosowania różnorodnych metod analizy danych i klasyfikacji, uzyskane rezultaty sugerują trudności w jednoznacznym modelowaniu i klasyfikacji danych dotyczących wartości odżywczych w menu Starbucks.